#include "FEMSpace_unstruct_dealii.h"
#include "Mesh.h"
#include "function2.h"
#include <ctime>

void generate_problem(EasyMesh*m, TemplateElement<2>*e, Equation<2>*equ,Domain<2>* domain  )
{
    Possion_2D_Dealii possionproblem(m, e, equ, domain);
    clock_t start, end;
    start = clock();
    possionproblem.setup_system();
    end = clock(); // 程序结束用时
    double endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "setup_system takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.assemble_system();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "assemble_system takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.dealwithBoundaryCondition();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "dealwithBoundaryCondition takes  " << endtime * 1000 << " ms"
         << endl; // ms为单位
    start = clock();
    possionproblem.solve();
    end = clock(); // 程序结束用时
    endtime = (double)(end - start) / CLOCKS_PER_SEC;
    cout << "solve takes  " << endtime * 1000 << " ms" << endl; // ms为单位
    int numdofs = m->getTotalNumDofs();
    dealii::Vector<double> Truesolution = dealii::Vector<double>(numdofs);
    vector<Dofs<2>> dofsList = m->getAllDofs();
    for (auto &dof : dofsList) {
        Point<2> pnt = dof.getPoint();
        double point[2] = {pnt[0], pnt[1]};
        int index = dof.getGlobalIndex();
        Truesolution[index] = u(point);
    }

    const DealiiSpMat& A = possionproblem.A();
    const DealiiSolution& solution = possionproblem.getSolution();
    const DealiiRHS& rhs = possionproblem.Rhs();
    dealii::Vector<double> residual(numdofs);
    double res_l2 = A.residual(residual, solution, rhs);
    dealii::Vector<double> err = solution;
    err -= Truesolution;
    cout << "relative error is: " << err.l2_norm()/Truesolution.l2_norm() << endl;
    cout << "the l2 norm of residual is:" << res_l2 << endl;
    // cout << "also l2 norm:" << residual.l2_norm() << endl;
    cout << "norm square:" << residual.norm_sqr() << endl;
    cout << endl;
};

void test_porder(int order) 
{
   vector<Point<2>> p = {{-1.0, -1.0}, {1.0, -1.0}, {1.0, 1.0}, {-1.0, 1.0}};
    Domain<2> *domain = new Rectangle2D(p);
    TemplateElement<2> *e;
    assert(order >= 1 && order <= 4);
    if (order == 1) {
        e = new Triangular_1_Element();
    } else if (order == 2) {
        e = new Triangular_2_Element();
    } else if (order == 3) {
        e = new Triangular_3_Element();
    } else if (order == 4) {
        e = new Triangular_4_Element();
    }
    Equation<2> *equ = new PossionEquation<2>();
    equ->SetRightHandsTermFunction(fe);
    equ->SetBoundaryConditionFunction({g1, g2, g3, g4});
    EasyMesh *m = new EasyMesh("D", order);
    m->distributeDofs();
    generate_problem(m, e, equ, domain);
    for(int i = 0; i < 1; i++){
      m->refine();
      m->distributeDofs();
      generate_problem(m, e, equ, domain);
    }
};

int main(int argc, char *argv[]) {
  int order;
  // cout << "Please input the order of Mesh:\n";
  // cin >> order;
  test_porder(1);
}
